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ABSTRACT 

In a cool neutron star (T < 10 6 K) endowed with a rather high magnetic field (B > 10 13 G), a phase 
transition may occur in the outermost layers. As a consequence the neutron star becomes "bare" , i.e. 
no gaseous atmosphere sits on the top of the crust. The surface of a cooling, bare neutron star does 
CO , not necessarily emit a blackbody spectrum because the emissivity is strongly suppressed at energies 

below the electron plasma frequency, oj p . Since uj p rs 1 keV under the conditions typical of the dense 
electron gas in the condensate, the emission from aT~ 100 eV bare neutron star will be substantially 
depressed with respect to that of a perfect Planckian radiator at most energies. Here we present a 
^ ' detailed analysis of the emission properties of a bare neutron star. In particular, we derive the surface 

emissivity for a Fe composition in a range of magnetic fields and temperatures representative of cooling 
isolated neutron stars, like RX J1856. 5-3754. We find that the emitted spectrum is strongly dependent on 
the electron conductivity in the solid surface layers. In the cold electron gas approximation (no electron- 
lattice interactions), the spectrum turns out to be a featureless depressed blackbody in the 0.1-2 keV 
band with a steeper low-energy distribution. When damping effects due to collisions between electrons 
and the ion lattice (mainly due to electron-phonon interactions) are accounted for, the spectrum is more 
depressed at low energies and spectral features may be present, depending on the magnetic field strength. 
\^ ' Details of the emitted spectrum are found, however, to be strongly dependent of the assumed treatment 

£Sj ! of the transition from the external vacuum to the metallic surface. The implications of out results to 

CO ■ RX J1856. 5-3754 and other isolated neutron stars are discussed. 

Subject headings: radiative transfer — stars: individual (RX J0720. 4-3125, RX J1856. 5-3754) — stars: 
neutron — X-rays: stars 

CO . 

O 1 !• INTRODUCTION 

■ More than 20 X-ray sources currently associated with isolated neutron stars (INSs) show evidence for a thermal compo- 
Q-T nent in their spectrum, in many cases superimposed on a power-law high-energy tail. In the commonly accepted picture 
q the hard tail is produced by non-thermal processes in the stellar magnetosphere. On the contrary, the thermal component 
5^ ■ originates at the surface, while the star cools down and internal energy is progressively radiated away. Since these objects 
^ \ are not complicated by strong accretion signatures, detailed observations of their thermal component provide a powerful 
C3 • tool to investigate directly the properties of the neutron star (NS). If the thermal emission originates in the NS atmo- 
L| ' sphere, the detection or absence of spectral lines and edges may constrain the chemical composition and/or magnetic field 

. t-H ! through a comparison with computed models (see e.g. Shibanov et al. 1992; Rajagopal, Romani & Miller 1997; Pons et 

■ al. 2002). Furthermore, empirical insights may be derived for the NSs mass, radius and equation of state (e.g. Lattimer 
Cj ; & Prakash 2001). 

The family of thermally emitting INSs includes seven peculiar objects serendipitously discovered in ROSAT PSPC 
pointings (see e.g. Treves et al. 2000 and Motch 2001 for reviews; Zampieri et al. 2001). These sources (hereafter 
referred to as ROSAT INSs) are characterized by remarkably similar properties among which are: soft, thermal spectrum 
with kT ~ 100 eV; low X-ray luminosity, Lx ~ 10 30 — 10 31 ergs _1 , and low column density, Njj ~ 10 20 cm~ 2 ; no 
association with a supernova remnant; pulsations in the 5-20 s range (detected in four sources so far). Until recently 
the spectral properties of the seven ROSAT INSs were not known in detail. PSPC observations provided evidence that 
a blackbody spectrum gives a satisfactory description of the data in all cases. However, the scanty statistics prevented 
definite conclusions being drawn about the X-ray spectral energy distribution (SED), and did not allow for the detection 
of spectral features. While this situation has not yet improved for the fainter sources, the two brightest ROSAT INSs, 
RX J1856. 5-3754 and RX J0720. 4-3125, have been the target of deep observations with Chandra and XMM-Newton. 
RX J0720. 4-3125 was observed with XMM-Newton for the first time on May 2000, in a 62.5 ks pointing (Paerels et al. 
2001; Cropper et al. 2001). The EPIC-PN spectrum is well-represented by a blackbody and no spectral features have been 
detected, apart from variations in the column density with pulse phase that may be explained in terms of energy-dependent 
beaming effects or cyclotron absorption (Paerels et al. 2001; Cropper et al. 2001; Haberl et al. 2003). The X-ray flux 
shows a modulation with a period of 8.31 s and a pulsed fraction of ~ 15% (Haberl et al. 1997; Cropper et al. 2001). Using 
ROSAT and XMM-Newton data, Zane et al. (2002) were able to derive the period derivative (P~5x 10 _14 ss _1 ) which, 
when interpreted in terms of magneto-dipolar braking, implies a surface magnetic field of ~ 2 x 10 13 G. The case of RX 
J1856. 5-3754 is even more striking. A 50 ks Chandra LETGS observation has convincingly shown that RX J1856. 5-3754 
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appears to have a featureless X-ray continuum, for which a simple blackbody yields a better fit than more sophisticated 
atmosphere models (Burwitz at al. 2001). Analysis of a very long (~ 500 ks) Chandra observation obtained in October 
2001 further reinforces this conclusion and rules out the presence of strong spectral features (Drake et al. 2002). However, 
the presence of broadband departures from a pure blackbody spectrum have been claimed by Burwitz et al. (2003) in 
both XMM-Newton and Chandra data. No X-ray pulsations have been detected and the upper limit on the pulsed fraction 
is now down to < 1.3% (Ransom, Gaensler & Slane 2002; Drake et al. 2002; Burwitz et al. 2003). RX J1856.5-3754 was 
the first source in this class for which an optical counterpart has been found (Walter & Matthews 1997; van Kerkwijk 
& Kulkarni 2001a). When used in conjunction with Chandra data, the recently measured star distance (about 120 pc; 
Kaplan, van Kerkwijk & Anderson 2001; Walter & Lattimer 2002) yields a radiation radius of only <~ 5-6 km (Drake et al. 

2002) . Very recently Haberl et al. (2003) reported the discovery of spectral feature(s) in the EPIC/RGS data of RBS 1223. 
The feature, seemingly present also in Chandra data, is a deficit of counts with respect to the best-fitting blackbody at an 
energy of ~ 0.2 — 0.3 keV, and is quite broad. Its nature is still uncertain but it could be a proton cyclotron absorption 
line in a magnetic field ~4x 10 13 G, as discussed by Haberl et al. (2003). 

The small apparent radius and the blackbody X-ray spectrum led to the intriguing suggestion that RX J1856. 5-3754 
might be a strange quark star (Drake et al. 2002; Xu 2002). One of the motivations for such a claim is that bare quark 
stars, i.e. those not covered by a layer of hadronic matter, do not have an atmosphere and would presumably emit a 
pure blackbody spectrum (as suggested for instance by Xu 2002). However, further investigations are definitely needed to 
assess the spectral properties of these objects. While a quark star may be a conceivable option (see however Thoma et al. 

2003) , present observations of RX J1856. 5-3754 do not necessarily demand this solution and more conventional scenarios 
involving a neutron star are certainly possible. Neutron star models based on a two-temperature surface distribution can 
account for both the X-ray and optical emission of RX J1856. 5-3754, giving at the same time acceptable values for the 
stellar radius (Pons et al. 2002; Walter & Lattimer 2002; Braje & Romani 2002). However, the problem of producing 
a featureless spectrum from a NS has not been consistently solved yet, although possible ways to suppress the spectral 
features which are expected in optically thick atmospheric models have been outlined (see e.g. Braje & Romani 2002). 

In this paper, we consider an alternative explanation for the peculiar X-ray spectrum of RX J1856. 5-3754. In the 70's 
it was commonly accepted that radiation emitted by NSs came directly from their solid surface (e.g. Brinkmann 1980 and 
references therein). Later, the role of the thin gaseous layer which covers the star crust in shaping the emergent radiation 
spectrum was appreciated and model atmospheres became the standard tool for interpreting the observed emission from 
isolated NSs. However, highly magnetized neutron stars may be left without an atmosphere if they are cool enough. 
The reason for this is the onset of a phase transition that turns the gaseous atmosphere into a solid when the surface 
temperature drops below a critical temperature T cr i t , which, for a given chemical composition, depends on the stellar 
magnetic field (see Lai & Salpeter 1997 and Lai 2001 for a recent review). The determination of T cr i t is still uncertain and, 
in particular, only preliminary calculations are presently available for heavy element (such as Fe) surface compositions 
(Lai 2001 and references therein). In §2 we show that, given the large uncertainties on the conditions for Fe condensation, 
it is possible that RX J1856. 5-3754 and (marginally) RX J0720. 4-3125 have surface temperatures below the critical value 
and may be then "naked" or "bare" neutron stars. 

The idea that RX J1856. 5-3754 might be a solid surface NS was suggested earlier by Burwitz at al. (2001) and Zane 
et al. (2003) (see also Burwitz et al. 2003). As first discussed by Lenzen & Trumpcr (1978), a severe reduction in the 
NS surface emissivity occurs at energies below the electron plasma frequency. Under the conditions for which a ROSAT 
INS is bare (T < 100 eV, B > 10 13 G), the plasma frequency in the surface layers corresponds to energies > 1 keV 
(see § 2), so the NS is expected to radiate less efficiently than a blackbody emitter at soft X-ray energies and below. 
This is of great potential importance, since it might help reconcile the observed radiation radius with current theoretical 
predictions of NSs radii. Also, it may hold the key for explaining the featureless blackbody spectrum observed in some 
of these sources. In § 3 we address in detail the question of the nature of the surface emissivity of a bare NS with a pure 
Fe composition. The method we use is similar to that employed by Brinkmann (1980), who was the first to investigate 
this issue in connection with X-ray pulsars. We found that the emissivity, and hence the shape of the emitted spectrum, 
depends crucially on the conductivity of the star crust. In § 3.1 we analyze a simple (albeit unrealistic) model in which 
only the contribution of a cold electron plasma to the dielectric tensor in the star interior is accounted for. Results for this 
case are qualitatively similar to those of Brinkmann (1980). Proper account for the damping produced by interactions 
of (degenerate) electrons with the ion lattice (mainly through electron-phonon collisions, e.g. Potckhin 1999), however, 
introduces qualitative changes to the above picture, as is discussed in § 3.2. The relevance of our model to RX J1856. 5-3754 
is finally discussed in § 4. 

2. BARE NEUTRON STARS 

In this section, we explore the possibility that some of the cooler isolated neutron stars (T < 100 cV) are left without an 
atmosphere by a phase transition in the surface layers at large magnetic fields (B > 10 13 G). We point out that, although 
this is unlikely for a light element (H, He) composition, it might be the case for heavy elements (such as Fe), at least 
for some sources notably including RX J1856. 5-3754. If indeed some INSs are bare, the question of the nature of their 
emitted spectrum arises. In a neutron star with metallic surface layers (here and in the following Z and A denote the 
atomic number and weight of the constituent element) the density at zero pressure is given by (e.g. Lai 2001) 



Ps » 560 AZ-V 5 B% 5 gem" 3 , 
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where B\2 = B/10 G. The electron plasma frequency is then 



hLUp = h J^i „ 0.7^/^3/5 f p_ \ 1/2 keV; (2) 




where n e is the electron density. In the following we take as a reference value for the plasma frequency that given by 
eq. (2) with p = p s and Z = 26; w Pi o represents then the plasma frequency in a pure iron medium with density p s assuming 
that all electrons are in the conduction zone (see §3.1 for a further discussion). 

Cool NSs (T < 100 eV) emit most of their thermal radiation below the plasma frequency and substantial deviations 
from a pure blackbody spectrum are expected as a result of the large absorption at uj < v p . Spectral features should also 
appear around the electron cyclotron frequency at 

eB 

Hlu b = — ~ 11.6 B 12 keV , (3) 

m e c 

but, since we focus on field values B > 10 13 G, these fall well outside the X-ray range accessible to the Chandra LETGS 
and XMM-Newton EPIC-PN, and are of no immediate interest. 

The properties of atoms and condensed matter are qualitatively changed by magnetic effects when b = B/B 3> 1, 
Bo = m e e 3 c/h 3 ~ 2.35 x 10 9 G. Theoretical research on matter in superstrong fields started over 40 years ago and 
although many uncertainties still remain, much progress has been made, especially for H and He compositions (see Lai 
2001 and references therein). In particular, when 6 > 1, electrons are strongly confined in the direction perpendicular 
to the magnetic field and atoms attain a cylindrical shape. Moreover, it is possible for these elongated atoms to form 
molecular chains by covalent bonding along the field direction. Interactions between the linear chains can then lead to 
the formation of three-dimensional condensates. As discussed by Lai & Salpeter (1997) and Lai (2001), in the case of 
hydrogen the infinite linear chains (and metallic hydrogen) are certainly bound, favoring the possibility of condensation 
for sufficiently low temperatures and/or strong magnetic fields. The critical temperature below which phase separation 
between condensed H and vapor occurs is 

T c H rit « O.lQoo (4) 

with 

Qoo w 194.1B? 2 37 - 4.4 (InB 12 - 6.05) 2 - tkv p , p - J (u% tP + wjjj 1/2 + \hu) B , v eV . (5) 

where w PiP and uj b _ p are the proton plasma and cyclotron frequencies. 

For heavier elements (such as Fe), the lattice structure and the cohesive properties of the condensed state are very 
uncertain and are different from those of H and He. For instance, unless the field is extremely high (-B12 3> 100), it 
is likely that the linear chains are unbound for Z > 6. More recent computations of the cohesive energy Q s of the 3D 
condensate showed that Q s is only a tiny fraction (~ 0.5%) of the atomic binding energy, correcting earlier overestimates 
(Jones 1986, see also Neuhauser, Koonin & Langanke 1987): 

Q s <0.05|Kw|~^ 9/5 -B 2 2 /5 eV forZ>10. (6) 

On the other hand (see again Lai 2001), even such a weak cohesion of the Fe condensate can give rise to a phase transition 
at sufficiently low T. The critical temperature at which phase separation occurs can be estimated by equating the ion 
density of the condensed phase near zero pressure [eq. (1)] to the gas density in the vapor (Lai 2001) 

Pg » 390A 5 / 2 T 5 / 2 exp (-Q s /T) g cm" 3 . (7) 

This gives 

T^ t <0.1Q s ^27B% 5 eV. (8) 

It should be noted that, although representing the more recent available estimates, these expressions for heavier elements 
are still quite crude: all models are approximate near zero pressure and the structure itself of the lattice is very uncertain. 
For our purposes, they should be regarded as being typically accurate to an order of magnitude. Also, the vapor density 
becomes much less than the condensation density and a phase transition is unavoidable only when the temperature drops 
below ~ T crit /2 (see Lai 2001). 

The critical condensation temperatures for H and Fe are plotted as a function of B in Figure 1. The filled circles 
show the position in the B-T plane of the coolest (T < 100 eV), thermally emitting INSs for which an estimate of the 
magnetic field is available (see table 1). We have also included RX J1856. 5-3754 in Figure 1; its position is indicated by 
a horizontal line since its magnetic field is not presently known. In order to obtain the local surface temperature, i.e. the 
quantity reported in Figure 1, a gravitational red-shift correction was applied to the values listed in table 1, according to 
the expression T sur f = (1 + z)T\,b where (1 + z) _1 = yjl — 2GM/c 2 R ~ 0.8 (M, R are the star mass and radius). Here 
Tbb is the color temperature, as derived from the blackbody fit. 

It is apparent from Figure 1 that all INSs have a temperature well in excess of the H critical temperature: if surface 
layers are H-dominated, the presence of a gaseous atmosphere is unavoidable. On the other hand, if INSs have not 
accreted much gas, one might expect to detect thermal emission directly from the iron surface layers. If this is the case, 
the outermost layers of RX J1856. 5-3754 (depending on the magnetic field), and possibly RX J0720. 4-3125, might be in 
form of hot condensed matter, in which case the usual radiative transfer computations do not apply. 
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3. THE SURFACE EMISSIVITY 

The emission properties of the neutron stars surface have been first analyzed by Lenzen & Triimper (1978) and in some 
more detail by Brinkmann (1980, hereafter B80). Both of these works were aimed to X-ray pulsars, where the surface 
temperature is a few keV's, and treated the medium inside the star as a cold electron plasma, neglecting all possible effects 
due to electron degeneracy and ion lattice (primarily through electron-phonon interactions). The (constant) damping 
frequency which appears in B80 calculations is mainly used to smear the resonance at the cyclotron frequency. Moreover, 
birefringence in the magnetized vacuum outside the star was not accounted for. In this section we derive the NS surface 
emissivity following an approach similar to that discussed in B80. To better illustrate the importance of electron-phonon 
interactions, we first consider a pure cold electron plasma, repeating Brinkmann's calculation for the parameter values 
appropriate to cold isolated NSs (§3.1). A complete treatment which includes the polarization properties of magnetized 
vacuum is presented in Appendix A. Since, as we show there, this more general approach is quite cumbersome and only 
gives tiny differences with respect to the simpler one based on unpolarized radiation, the latter is used below. In §3.2 we 
analyze the more realistic case in which the damping of electromagnetic waves produced by the presence of the ion lattice 
is included. 



3.1. The Cold Plasma Case 

We start considering the medium inside the star as a cold electron plasma and neglect the damping of free electrons 
due to collisions. We introduce a cartesian frame as in B80 (see his figures 1 and 2) with the z-axis parallel to the surface 
normal. The direction of the incident wave vector k is specified by the angle of incidence i and the azimuth (3. The 
magnetic field direction b = B/B is at an angle a with respect to the z-axis and b lies in the x — z plane. Given a star 
surface element dA — 2-kR 2 sin 9 d8 at magnetic co-latitude 9, we first compute the total reflectivity p^ of the surface for 
incident unpolarized radiation. Then, since the absorption coefficient is simply — 1 — p u , Kirchhoff's law yields the 
emissivity j u — a u B u (T), where T is the temperature of the emitting element. In general, p u depends on the direction 
of the refracted ray (see below). Therefore, the monochromatic flux f u emitted by the surface element must be computed 
by integrating over all incident directions, 

U= / j ul (i,0,6)smidid(3. (9) 

JO J-ir/2 

The flux emitted by the entire surface is given by 4 

i r i f-27t />tt/2 

F " = T152 UdA=- sm6d6 / j u (i, 0, 9) sim di d(3 . (10) 

47TK J sphere Z J Jo J-tt/2 

At the surface, an incident electromagnetic wave, described by its electric field E and wave vector k, is partly reflected 
(E",k") and partly refracted. Due to the birefringence of the medium, the refracted wave is the sum of an ordinary 
(E' 1; k[) and an extraordinary (E^, k' 2 ) mode. In order to compute the reflectivity, we need to solve the dispersion relation 
and to compute the refractive index n for the two modes of propagation. In our frame, the dielectric tensor for a cold 
electron plasma is given by 

S cos 2 a + P sin 2 a — iD cos a sin a cos a(P — S) 

iDcosa S —iDsina j (11) 

sinacosa(P — S) iDsina P cos 2 a + S sin 2 a 



with 



,2 



P=l-^f, (13) 
(R±L)/2. (14) 



S 
D 

By introducing the Maxwell tensor Ay = k^kj — k' 2 5ij + {uo 2 /c 2 )eij , where k[ are the cartesian components of k' and 
k 12 = k^kl, the dispersion relation is obtained by imposing |Ay| = 0. For our purposes it is convenient to write the resulting 
expression in terms of angle of incidence i, and the (complex) refractive index n = k'c/u>. By using an expression formally 
analogous to Sncll's law n — sin ij sin <d (where now O is a complex quantity which replaces the angle of refraction while 
i is real; see e.g. Marion 1965), it is 



n 4 (P + v sin 2 a) + n 2 (gv - 2PS + u sin 2 a) + PRL + gu = 
sinisin(2a) cos/3(n 2 — sin 2 z) 1 / 2 (u + n 2 v) ; 

4 Viewing angle effects have been neglected in evaluating the flux from eq. (10). 



(15) 
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in the previous expression v — S — P, u = PS — RL and g = sin 2 i[l — sin 2 a(l + cos 2 (3)}. Squaring eq. (15) gives a 
fourth order polynomial equation in n 2 which can be solved analytically. Clearly only two out of four solutions satisfy the 
original dispersion relation and represent the refractive indices for the two propagation modes in the magnetized plasma, 
n m , m = 1 , 2. As noted by B80, the only practical way of finding the two meaningful roots is to substitute them back 
into eq. (15) and check numerically the residual. This, however, turned out to be troublesome for some values of the 
parameters, as we discuss later on. For i = 0, a = or ir/2, = tt/2 or 3tt/2, the right hand side of eq. (15) vanishes and 
the dispersion relation reduces to a quadratic equation in n 2 which is then solved instead of the quartic. 

Once the refractive indices are known, we can solve the wave equation for the two refracted modes Xij(n m )E' m j = 0, 
where E' m - are the cartesian components of E' m , obtaining the two ratios E' X /E' and E' /E' z . We performed the 
calculation (double-checked with the aid of an algebraic manipulator) , obtaining 



E 1 

m.x 



EL 



= a r , 



-r? m sin 2 i sin cos (3 — iD sin 2 i cos a + iD cos f3 sin a sin i 



riL — sin 2 i 



— sin/3 (P — S) sin a cos a sin i \J — sin 2 i 

+ sin 2 i sin j3 cos (3 (P cos 2 a + S sin 2 a) + iD cos aP] 

x — ri 1 m sin i \J n 2 ^ — sin 2 i sin (3 + iD sin ar^ m — iD sin a sin 2 i cos 2 (3 

— iD cos a cos 0s\x\i\J — sin 2 i + sin i \J — sin 2 i [sin (3S 



+ sin (3 sin 2 a (P — S)] — (P — S) sin a cos a sin 2 i sin (3 cos (3 — iD sin aP] 



(16) 



E' 

m,y 

~E' 



(sin 2 i sin (3 cos [3 — iD cos a) + sin (3 sin i \J n 2 ^ — sin 2 i + iD si 



sin a 



x (sin 2 (3 sin 2 i - n 2 m + S) 

While the previous expression for b m agrees with that given in B80, our result for a m is different and we were unable to 
recover his expression. 

The ratios in eq. (16) are then inserted into the Fresnel equations which fix the boundary conditions at the interface 
between the two media (see Jackson 1975 and cqs. [17]-[18] in B80). This allows the derivation of the components of the 
electric field of the reflected wave parallel and perpendicular to the plane of incidence (E^ and E'{_) in terms of the same 
components of the incident wave {E\\ and E±). We re-derived the expressions given in B80 and found 



E = V~ 



E'l=V~ 



B X B 2 (1 

1 — W\ 
1 + W\ 



f w 2 ) (1 + wi 
1 + w 2 



■E, 



+ (A_-B_)E ll 
E± + 2BiB 2 {wi 



Wi) E\\ 



(17) 



where 



A± =B 1 (l + w 1 ) 
B ± =B 2 (l + w 2 ) 



A 2 ( w 2 cos i ± 
Ai I w\ cos i ± 



COSi 

1 

COSi 



sm i 



(18) 



' n — sin 2 i I cos i , A m = b m sin (3 — a m cos 0, B m = b m cos + a m sin and T> = A + — B + . These expressions 

were double-checked with the aid of an algebraic manipulator and differ again from those in B80; we also note that our 
definition of B± is different from that of B80. 

The reflection coefficient for unpolarized radiation can be expressed as the combination of the reflectivity of parallel 
and perpendicularly polarized incident waves, p w = {p\\ tU + p_L,v)/2. Since the reflectivity is defined as the ratio of the 
reflected to the incident wave amplitudes, pi\ u is the sum of the square moduli of the coefficients of En in eqs. (17). 
Similarly, p± tLU is obtained by adding together the square moduli of the coefficients of E±. 

The absorption coefficient a u = 1 — has been computed numerically in the relevant angular ranges following the 
procedure outlined above and the results have been used to evaluate the integral in eq. (9). Although the numerical 
scheme is rather straightforward, care should be used since the refractive index becomes resonant where the coefficient of 
the higher order term in eq. (15) vanishes (B80; see also Melrose 1986). This happens at P + vsin 2 a = 0, that is to say 
at the two frequencies 



■CJf 



1 ± 



1 



cos a 



,2 \2 



1/2' 



(19) 



Assuming no collisional damping has the main advantage that the fourth order polynomial obtained by squaring 
the dispersion relation has real coefficients and its roots are either real or complex conjugates in pairs. Numerical 



6 



experimenting shows that the roots develop an imaginary part only close to the resonances and are real and distinct 
otherwise. As discussed by Melrose (1986), there exist two cut-off energies at which one of the two indices vanishes. The 
cut-off energies are a-dependent and the smallest one, u , falls in the range < uj n < uo + . This corresponds to the 
appearance of an evanescent mode, which can not propagate into the medium. Modes for which the refractive index has 
a large imaginary part are severely damped (e.g. Jackson 1975), so they can not penetrate much below the surface either. 
We point out that the existence of these damped waves is not in contradiction with having neglected collisional damping. 
The "conductivity" afj of a cold plasma can be computed from the dielectric tensor (11) using the standard relation 

47T i 

£ij =S ij+ i— o% (20) 

and is therefore purely imaginary (see Jackson 1975; Meszaros 1992). Quotation marks are placed on "conductivity" 
because there is no resistive loss of energy in this case. Yet, depending on their frequency and grazing angle, electromagnetic 
waves may be exponentially damped in the cold electron plasma. 

In all cases in which one of the acceptable roots is real and negative or it is complex and its imaginary part exceeds a 
given value, only one refracted mode survives. Consequently, we adopt a "one mode" description to derive the reflectivity. 
By specializing the previous calculation to a single mode (labeled '1' for convenience), and defining C\ = £3101 + 63261+633, 
we obtain (see also Jackson 1975) 



d arctan -A V j ~ sin 2 i cos i 

E W d arctan, + A 11 ' E± I 2 = TT~ ~ E± (21) 

1 1 J nf - sin 1 + cos 1 

from which the reflectivity follows. 

As we mentioned earlier, selecting the two relevant roots of the fourth order polynomial may become difficult in some 
parameter ranges. While for most values of the energy and angles, the two meaningful roots produce a residual many 
orders of magnitude below that of the spurious ones, in some cases all residuals are small and of the same order. We 
encountered situations in which, owing to round-off errors and despite the use of quadruple-precision complex arithmetic, 
one of the relevant roots produced a residual larger than that of the spurious solutions. A bad choice of the roots usually 
gives negative absorption. In these cases the calculation is repeated with a different choice of the roots until is positive 
and its value close to those computed for neighboring values of the parameters. 

The quantity f^/B^iT) (sec cq. [9]) is shown in Figure 2 as a function of energy for B = 10 12 G and different values of the 
angle a. Here the surface temperature is assumed constant. In general, the emissivity becomes lower as the magnetic field is 
more and more tilted with respect to the surface normal. The strong absorption around the resonant frequency w_ is clearly 
seen; the absorption dip becomes more pronounced (and the surface emissivity decreases) as a approaches n/2. The total 
monochromatic emissivity is obtained integrating over the entire stellar surface (eq. [10]), once the magnetic field topology 
is specified. Here and in the following we assume that the field is dipolar, B = B p [(4 — f) cos 2 9 + /] 1 ^ 2 /2, where B p is the 
polar field strength and / ~ 1.2 accounts for general-relativistic corrections in a Schwarzschild space-time (see e.g. Pavlov 
& Zavlin 2000). Accordingly, a is related to the magnetic co-latitude 9 by cos 2 a = {4 — 4//[(4 — /) cos 2 9 + /]}/(4 — /). 

The quantity F^/B^iT) is shown in Figure 3 for different values of B p (models computed accounting for e-phonon 
interactions are also shown, see § 3.2 for details). Again, the surface temperature is taken constant. As expected, 
integration over 9 smears out any strong feature around w_, as it can be seen comparing the full line in Figure 3 with 
those in Figure 2 (for B p > 10 12 G the angle-averaged resonant feature lies outside the energy range that we have 
considered). Below u>_ one of the two modes is non-propagating, a so-called whistler. Whistlers have a very large 
refractive index, which diverges at and for lo — ► 0, and this explains the high reflectivity at these frequencies (e.g. 
Melrose 1986). The increase of the plasma frequency with B is responsible for the lower emissivity at larger fields and, 
if we restrict to energies below the (angle-averaged) resonance w_, the dependence on energy is about the same (i.e. the 
curves are nearly self-similar, see again Figure 3). 

The model presented so far has been computed using for the plasma frequency in the surface layers our reference value 
iv P fl. It should be stressed, however, that this value, which relics on the zero-pressure surface density given by eq (1), 
is just an estimate (see e.g. Lai 2001). In order to assess the effects of this on the emitted spectrum, we explored the 
parameter space by varying the electron plasma frequency around uj p ,o; this accounts also for the uncertainties on the 
number of free electrons per nucleon, i.e. deviations from our reference value Z = 26. 

While the energy dependence of the absorption coefficient does not change significantly, an increase in the plasma 
frequency at fixed B produces an overall decrease of the emissivity. This is illustrated in Figure 4, where the ratio of the 
emitted to the blackbody power in the 0.1-2 kcV band is plotted against the plasma frequency for B p = 3 , 5 x 10 13 G. As 
can be seen from Figure 4, in this simple description we expect that the star surface radiates only ~ 30% of the blackbody 
power if the plasma frequency (the density) is about a factor 6 (36) higher than the estimate provided by eq. (2) (eq. 
[1]). Also, since in this simplified description f u /B u (T) is found to vary with the magnetic field angle (Figure 2), viewing 
effects may be relevant with even larger depression expected if the star is viewed equator-on. 

3.2. Electron- phonon Interaction 

In the calculation presented so far all collisional damping effects have been neglected, therefore the only source of 
"conductivity" is the cold plasma. This is of course an oversimplification. In the solid crust of a NS with B ~ 10 12 — 10 13 G 
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and T ~ 10 6 K electrons are strongly degenerate (T <C Tp, where Tp is the Fermi temperature). Furthermore, if p ~ p s 
quantum effects due to the magnetic field are not negligible. In the degenerate surface layers charge (and heat) is 
transported primarily by electrons. Several efforts have been devoted to an accurate determination of the electrical 
conductivity a, mostly because this is the basic quantity governing the magnetic field evolution (see e.g. Flowers & 
Itoh 1976; Yakovlev & Urpin 1980; Potckhin 1999 and references therein). In particular, it is well known that while 
at temperatures above the crystallization temperature of the ions the main factor governing the electrical conductivity 
is scattering off ions, below the crystal melting point the dominant process is scattering by crystal lattice vibrations 
(phonons) through Umklapp processes. Eventually, if the temperature decreases further, impurities in the crystal structure 
and scattering off lattice defects start to be important. 

Magnetic fields in the neutron star crust complicate electron transport, in particular, making it anisotropic (Kaminker 
& Yakovlev 1981; Yakovlev 1984; Hernquist 1984; Potckhin 1999). Electrons move freely parallel to B, but their motion 
perpendicular to the field is quantized. Therefore, transport properties will be affected through the influence of the 
discrete spectrum on the density of states and collision times. The electrical conductivity tensor is computed from the 
transport tensor and can be expressed as (Hernquist 1984; Potekhin 1999) 

f 2 ^b (e) . , ( df \ 

where e is the electron energy (including the rest energy m e c 2 ), / (e) is the Fermi-Dirac distribution, 

N B (e) = meLJB 2 V" g n (e/c) 2 - (m e c) 2 - 2m e huj B n , (23) 

g n is the statistical weight of the n'^-Landau level and n max is the maximum Landau number for a given energy e (see 
Potekhin 1999 for all details). The quantities (e) in eq. (22) play the role of relaxation times. The tensor has three 
independent components: one parallel to the field {a zz = an), one transverse (a xx — a yy = a±) and one off-diagonal 
(Hall) component a xy = cfr, which is non-dissipative. In general evaluating eq. (22) requires energy-integration, but for 
strongly degenerate electrons and not too close to the Landau thresholds this becomes unnecessary. In this case it is 
<7ij w (e 2 n e c 2 / 'ef) {sp) where ep is the Fermi energy. 

The way in which the magnetic field affects the charge transport depends on its strength. According to Potckhin et 
el. (2003), the non-quantizing (classical) regime occurs if T 3> Tb, where kT B = heB/(m e c). Even a non-quantizing 
magnetic field, which essentially does not affect the thermodynamic properties of matter, hampers transverse motion and 
produces Hall currents. In the opposite case, T -C Tb, one can distinguish between a weakly quantizing (when electrons 
populate several Landau levels) and a strongly quantizing regime (when the magnetic field confines most electrons in 
the ground Landau level). Transition to the strongly quantizing regime occurs below the first Landau threshold, i.e. for 

p < p B « 7 x 10 3 (A/Z)B% 2 g/cm 3 . Note that the zero pressure density p s is much smaller than p B for B ;§> 10 10 G. 

So far, the most complete expressions for the electrical conductivity are those computed by Potekhin (1999), by taking 
into account correlation effects in the strongly coupled Coulomb liquid and multi-phonon scattering in the Coulomb 
crystal (Baiko et al. 1998). These results show that in a weakly quantizing field the conductivity in the longitudinal and 
transverse direction oscillates around their classical values. However, such oscillations are quite prominent in the regime 
of strong quantization, where they may reach several orders of magnitudes. In this case the transport properties of the 
matter are very different from those in the classical regime. 

In order to include effects of electron-phonon scattering in our computation, it is convenient to introduce the effective 
relaxation times. In the regime of strongly degenerate electrons, when thermal averaging is unimportant, the relaxation 
times are related to the conductivity along the various directions by t xx = r yy = Ana^/ujp, r zz = Aira^/ijjp (e.g. Ziman 
1978; Yakovlev & Urpin 1980; Potekhin 1999). Then we proceed in a standard (although approximate) way by summing 
the relevant effective frequencies of the different processes. Note that, while for transport along the field the effective 
collision frequency for electron-phonon damping is simply ~ r,r = t^ 1 , for transport across the field it is not always 

~ t xx . In the strong field regime it becomes directly proportional to t xx . Physically this reflects the fact that the 
relaxation time in the transverse direction is, in the strong field regime, longer than the time between electron-phonon 
collisions (tj_) which in turn determines the damping according to Heisenberg's principle. Accordingly, we derived the 
collision frequency from the interpolation formula by Potckhin (1999), valid at any field strength t xx = tj_/(1 + ^b t i_ )• 
We focus on the case B || z (i.e. a = 0) and write the conductivity tensor of the pure plasma component (eq. [20]) in 
rotating coordinates (e+,e_,e z ), e± = ± ie y (Meszaros 1992). In this frame the conductivity tensor a pl is diagonal 5 

1 + 

- 1 , • W 
V 1/ 

5 The plasma polarization tensor, whose elements are directly related to the collision times, is in fact Hij = v^ 1 (Sij — tij) = —(Hn/vw)<Tij, 
where v = lo^/lo 2 . In a magnetic field and in rotating coordinates with B || 2, the polarization tensor is diagonal (see e.g. Meszaros 1992). 
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where 



The conductivity tensor related to electron-phonon scattering can be written (again in rotating coordinates) as 



~e— 



Wc then obtain 



-tot 



-tot 

yy 



where the quantities 



Up , 
u yy 




(25) 



(26) 
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47TW 1 + Ub/u + tLO^/u) 
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(28) 



play the role of damping frequencies in the two different directions. By transforming back to non-rotating coordinates we 
finally get 

L tot + R tot - 2 i{L tot -R tot ) 



i4ir 



-i (L tot - R tot ) L tot 




Jftot 





2 

<2ptot 



from which the dielectric tensor follows as 



s** ot 


_ iD tot 







gtot 











ptot 



(29) 



(30) 



In the previous expressions, the off-diagonal terms correspond to the (non-dissipative) Hall conductivity, and S tot , D tot 
are defined as 5, D in § 3.1 but with R, L, P replaced by 

,2 

(31) 



L tot 



= 1 _ UJ P 



Ul 



D 



UJ^ LU =F LOB + 
,2 



ptot _ ^ 



P 



w + iwjP 



(32) 



A further rotation by an angle a accounts for the misalignemcnt between B and z and gives the dielectric tensor e'° in 
the same form as in eq. (11). 

We have repeated the computation of the monochromatic absorption coefficients by following the same method as in 
§ 3.1, but using the dielectric tensor e*°*. The conductivites <J\\.±_ have been computed numerically for the appropriate 
values of B, p and T 6 . The inclusion of electron-phonon damping seriously affects the emission properties of the surface, 
as can be seen from Figures 3 and 5 where the angle-averaged and angle-dependent emissivity is shown. 

Below ~ 1 keV the emissivity declines quite rapidly with decreasing photon energy. There is now a strong dependence 
on the magnetic field strength, with the suppression being more pronounced at lower fields. This is mainly due to the role 



played by electron-phonon damping in the transverse plane which is quantified by the real parts of a 



tot fjtot 

u yy ■ 



These arc 



« LOpLO 1 ^ I {Aituj b ) and increase monotonically with decreasing B. In particular, while for B p ~5x 10 13 G the star surface 
can radiate down to a few tens of eV's, for B p < 5 x 10 12 G, a sharp edge appears at ~ 300 eV and the star surface 
behaves as a perfect reflector at energies below w 100 eV. An absorption feature close to the cutoff energy loq is now more 
clearly seen at w 300 eV for the field strengths reported here. Contrary to the behaviour of the resonant frequency w_, 
the value of the cutoff energy changes only very weakly with B (see Fig. 5). This implies that the feature survives even 

6 The package CONDUCT. FOR developed by A. Potckhin and available at www.iofTe.ru/astro/conduct/condmag.html has been used. 
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after integration over the entire surface has been performed (see Fig. 3). As in the cold electron gas considered in § 3.1, 
the emissivity is strongly dependent of the angle between the magnetic field and the surface normal (see again Fig. 5), so 
line-of-sight effects are expected to be important. 

The different behaviour produced by the inclusion of damping can be understood as follows. We start rewriting the 
dispersion relations (eq. [15]) in terms of the angle of refraction O, and then compare the refractive indices computed 
with and without the damping terms. We consider the case B || z (i.e. a = 0; see Melrose 1986). The ensuing equation 
is quadratic in n 2 and therefore simpler to solve. Results are shown in Figure 6 for 6 ~ 43° and two different magnetic 
field strengths. As can be seen from the bottom panel of Figure 6, when collisional damping is included, the r? m develop 
a substantial imaginary part over a wide range of energies which can easily extend up to an order of magnitude below 
The region of interest is larger for decreasing field strength and/or increasing grazing angle. This causes the angle- 
dependent absorption feature that is seen in Figure 5 at log E ~ —0.6, and large a. Moreover, below a certain energy 
threshold both indices develop a large imaginary part so that both propagation modes are substantially damped. Since 
the penetration depth is 5 = c/[oTm(n)], in our angle-dependent computation we neglect the contribution of modes with 
Im(n) = X/(2nS) > 0.01 (dashed line in Figure 6). 

Of course the most important consequence is that the results now depend in a crucial way on the choice of the rejection 
limit, which in turn is only one of the aspects of the uncertainties in the physics at the vacuum/surface interface. The 
relevance of varying the adopted limit for rejecting on the surface emissivity is shown in Figure 7. 

4. DISCUSSION 

This investigation has been motivated by recent X-ray observations of RX J1856. 5-3754 and RX J0720. 4-3125, the two 
brightest among the seven isolated neutron stars discovered by ROSAT. In particular, detailed Chandra observations have 
convincingly shown that RX J1856. 5-3754 has a featureless thermal spectrum for which a simple blackbody distribution 
seems to provide a better fit to X-ray data than more sophisticated atmospheric models (Burwitz at al. 2001; Drake et al. 
2002). RX J1856. 5-3754 has a firmly established optical counterpart (Walter & Matthews 1997; van Kerkwijk & Kulkarni 
2001a, similarly to RX J0720.4-3125 (Motch & Haberl 1998; Kulkarni & van Kerkwijk 1998). Accurate spectroscopy 
and photometry with combined Very Large Telescope (VLT) and Hubble Space Telescope (HST) data has shown that 
the UV-optical energy distribution closely follows a Rayleigh- Jeans tail (van Kerkwijk & Kulkarni 2001a). However, as 
originally noticed by Walter & Matthews (1997), the Rayleigh-Jeans tail of the X-ray best-fitting blackbody underpredicts 
the optical flux by about a factor 6 (Walter & Lattimer 2002), and this has been taken as suggestive of emission from 
regions on the star surface with different properties (see below). Deep HST observations have also revealed the presence 
of a bow-shock nebula in H Q around RX J1856. 5-3754 (van Kerkwijk & Kulkarni 2001b). Walter (2001), by means of 
HST observations, derived the proper motion and parallax of the star, and obtained a distance of about 60 pc. However, 
this distance was shown to be in error by recent re-analysis of the same HST data (Kaplan, van Kerkwijk & Anderson 
2001), and of these data augmented by further observations (Walter & Lattimer 2002). These studies place the source 
at about 120 pc, or twice the original distance. The simultaneous determination of the distance and X-ray flux (under 
the reasonable assumption that it comes from the stellar surface) makes RX J1856. 5-3754 unique in its class inasmuch it 
allows a direct estimate of the radiation radius 

Roo = 4.25 ( -^—] ( -^-) 2 km. (33) 
VIOOpcy V60eVy v ; 

Based on the location of RX J1856. 5-3754 in front of the R CrA molecular cloud, Drake et al. (2002) have shown that 
their derived neutral H column density of 1 x 10 20 cm" 2 limits the distance to < 170 pc. Taken at face value, the expression 
for the radiation radius above yields for this distance a value of at most ^8.2 km. Such a figure is incompatible with 
current bounds on the stellar radius (as measured by an observer at radial infinity) based on theoretical investigations of 
the equation of state of matter at ultra-high densities (EOSs; see, e.g., Lattimer & Prakash 2001), 12 km < < 17km. 
This discrepancy motivated the suggestion that RX J1856. 5-3754 might be a strange/quark star (Haensel 2001; Xu 2002; 
Drake et al. 2002; Gondck-Rosinska, Kluzniak & Stergioulas 2002). 

More conventional scenarios involving a neutron star have been discussed in connection with RX J1856. 5-3754. Pons 
et al. (2002) explored non-magnetized model atmospheres with different compositions (H, He, Fe, Si-ash) in order to 
reproduce the emission properties of RX J1856. 5-3754. H/He spectra are almost featureless but they deviate substantially 
from a blackbody, showing an excess at higher energies. As was already noted by Campana, Mereghetti & Sidoli (1997) 
and further remarked on by Pons et al. (2002), the fit with H/Hc atmospheric models yields a column density distance 
< 10 pc for standard values of the star radius, more than one order of magnitude below the parallax measurement. 

Pons ct al. (2002) considered both a uniform thermal distribution and a two-temperature surface model. For the 
adopted distance of 61 pc (Walter 2001), they concluded that the apparent radius is « 7-8 km for a uniform thermal 
distribution — too small for any EOSs. Owing to the reduced area of the hotter X-ray emitting region, two-component 
models provide larger values for R^ and may explain also the excess at optical wavelengths over the best-fitting X-ray 
blackbody. With a revised distance of 117 pc, Walter & Lattimer (2002) argued that both the single-component heavy- 
element (Fe and Si) and two-component blackbody models of Pons et al. (2002) yield acceptable values for the stellar 
radius. In particular, they claim that a two-temperature model in which a blackbody at T — 15 eV is emitted by a region 
with an angular diameter five times larger than the X-ray blackbody (T = 63 eV) can reproduce the multiwavelength 
SED. 
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These conclusions are not without problems. Although detailed spectral calculations were announced, the present results 
of Walter & Lattimcr (2002) rely on the assumption that the two regions on the stellar surface emit a pure blackbody 
spectrum. How to justify this assumption for a neutron star, however, remains unexplained. Moreover, the presence of a 
small, hot region on the star surface might be difficult to reconcile with the lack of pulsations in the X-ray flux, as stressed 
by Drake et al. (2002), especially in the light of the present very tight limits on the pulsed fraction of this source (< 1.3 %; 
Burwitz et al. 2003). All heavy element spectra calculated by Pons et al. (2002) exhibit a variety of emission/absorption 
features in the soft X-ray range. No evidence for such features is present in the Chandra and XMM-Newton data. 

More recently Braje & Romani (2002) readdressed the two-temperature surface model for RX J1856. 5-3754, pointing 
out that several effects (isothermality, magnetic smearing, rotation) may act in suppressing the spectral features from 
an extended atmosphere with heavy elements. In particular, they discuss in detail the role of rotation, showing that 
phase-dependent Doppler shifts in a rapidly rotating neutron star (P w 1 ms) wash out all features, leaving a nearly 
planckian spectrum. Although such a short period can not be excluded on the basis of present data, the detected periods 
of other thermally emitting INSs are in the range w 0.1-10 s, about two orders of magnitude larger. In the Braje & 
Romani (2002) picture the X-ray emitting region is kept warm by external heating, and the genuine surface temperature 
should correspond to the cooler blackbody at T <~ 15 cV. A millisecond period appears hardly compatible with the star's 
age of, as implied by conventional cooling curves, « 10 6 yr. Furthermore, the energetics of the bow-shock nebula implies 

1 /2 

P = 4.6 (-B/10 8 G) ms. A ms spin period therefore necessarily demands for a very low field star. Such a low field seems 
hard to reconcile with the limit on the age derived again from the bow shock energetics, (B/10 12 G)(r/10 6 yr) ~ 3 — 4. 

Although two-temperature models appear promising in explaining the multiwavelength SED of RX J1856. 5-3754, no 
conclusive evidence has been provided yet that a near blackbody, featureless spectrum can be emitted by an extended 
atmosphere covering the stellar crust. An alternative possibility, originally suggested by Burwitz at al. (2001) (see also 
Burwitz et al. 2003) and further explored here, is that RX J1856. 5-3754 may be a solid-surface NS. If this is the case, 
a severe reduction in the surface emissivity has to be expected at energies below the plasma frequency, according to the 
analysis of Lenzen & Triimper (1978) and Brinkmann (1980). The bare NS model may ease the radiation radius problem. 
In fact, denoting with /e the ratio of emitted to blackbody power and assuming emission from the entire star surface, 
the value of now contains an additional f E 1 factor with respect to that given by eq. (33). As expected, the reduced 
surface emissivity acts precisely in the same way as a reduced emitting area, requiring a larger star radius. In order 
to represent a viable option, the bare NS picture must conform to three basic requirements: i) the conditions for the 
appearance of a solid phase should be met, at least within the present uncertainties; ii) the X-ray spectrum emitted by 
the surface should be very close to blackbody in the 0.1-2 keV range to match Chandra/ XMM observations; and hi) quite 
low values of /e (~ 0.1) should be possible for i?oo to be in the range allowed by current EOSs. 

In the absence of a measured period and period derivative, the magnetic held of RX J1856. 5-3754 is still a mystery. 
Given a surface temperature of ~ 70 eV, the magnetic field of RX J1856. 5-3754 should be in excess of 10 13 G, more 
probably at least 3-5 x 10 13 G, for its surface layers of to be in the form of condensed iron (see § 2 and Figure 1). Although 
rather high, such a field strength is well below the magnctar range and is noticeably shared by another ROSAT isolated 
neutron star: RX J0720. 4-3125 (Zane et al. 2002). So, although no definite conclusion can be drawn, the possibility that 
RX J1856. 5-3754 is a solid-surface NS is real. 

In the light of the results presented in § 3, the remaining two points are much more of an issue. When computed 
accounting only for the cold electron gas (§ 3.1), spectra show indeed only small departures from a blackbody in the 
0.1-2 keV band. A typical example is shown in Figure 8, where the computed spectra are plotted together with the 
best-fitting blackbody (again in the 0.1-2 keV range) for two different temperature distributions on the star surface. 
Deviations from a blackbody are below 15-20%, and this would make possible to reasonably fit Chandra data. The value 
of Je depends on the magnetic field and on the surface density (see Figure 4). In order to reach Je ~ 0.3 (which produces 
an increase in radius of ~ 2) one must invoke a density « 50 times larger than the zero pressure value given by eq. (1). 
Although the latter is only an approximation, such a large departure might not be realistic. For p ~ p Sl the predicted 
increase in radius is ~ 15%, which is insufficient to reconcile the radiation radius with canonical theoretical predictions, 
at least when viewing angle effects are neglected (see also Thoma et al. 2003). 

These models are, however, unrealistic. The cold electron gas assumption is a poor approximation at low energies where 
damping by electron-lattice interactions becomes progressively more important. When collisional damping (computed 
following Potekhin 1999) is accounted for, we found that the surface emissivity is substantially depressed below m 1 keV 
with respect to the cold electron gas case. At low fields (< 5 x 10 12 G) virtually no emission is expected at energies below 
« 0.5 keV, while the decline at low energies is not so sharp for B > 10 13 G (see Figure 3 and the discussion at the end of 
§3.2). The emerging spectrum is shown in Figure 9 for two representative values of the polar magnetic field. Despite the 
spectra deviating quite strongly from a blackbody distribution at low energies, the fit with a blackbody in the 0.1-2 keV 
range is still acceptable, with maximal deviations typically below 20%. At the lower field strength shown in Figure 9 
(B p = 2 x 10 13 G), Je ~ 0.35 which would imply a radius larger than the pure blackbody radius by about a factor two. 
This is definitely larger than what is predicted by the cold electron gas models with p = p s , and may be enough to provide 
an acceptable value of the stellar radius. However, at least for the uniform temperature distribution, for such values of the 
polar field the absorption feature around w ~ 300 eV is clearly present in the spectrum (see again Figure 9). The feature 
is not so pronounced at larger fields but Je becomes higher (~ 0.45) making the radius a problem again. One has also to 
bear in mind that the spectra shown in Figure 9 have been computed for a fixed rejection threshold. As Figure 7 shows, 
the choice of this parameter (even within a factor of a few) has a crucial influence on the shape of emitted spectrum. 
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Apart from the considerable uncertainties in current modeling of the physics governing the phase transition (see §2 and 
Lai 2001 for a more detailed discussion), we remind the reader that our spectra have been computed under a number of 
simplifying assumptions. A thorough discussion of the limitations of this kind of approach can be found in Brinkmann 
(1980). The greatest uncertainties arise because of the assumption of a sharp transition from vacuum to a smooth metallic 
surface, neglecting the effects of the macroscopic surface structure. We assumed the surface is made of pure iron, but 
different chemical compositions, or the presence of impurities in the iron surface, may change the results. Inside the star 
we neglect the role of bound electrons and further effects produced by the dissipation of those waves that are rapidly 
attenuated within a skin penetration depth. 

Finally, from a different perspective and regarding the possible application of the present work, we point out that the 
calculation of the complex refractive indices below the plasma frequency presented here may substantially contribute to 
the determination of the photon thermal conductivities of ultramagnetized neutron stars. These are, in turn, important 
for the accurate calculation of the thermal structure and cooling of these objects (see Potekhin et el. 2003). 

We are deeply indebted to A. Potekhin for a critical reading of the manuscript, for his comments and suggestions and 
for pointing out the relevance of this work in connection with photon thermal conductivities computations. We thank 
M. Chieregato and K. Wu for many helpful discussions during the earlier stages of this work. We also acknowledge 
an anonymous referee whose penetrating questions and comments greatly improved an earlier version of this paper. 
This work was partially supported by the Italian Ministry for Education, University and Research (MIUR) under grant 
COFIN-2002-027145. JJD was supported by NASA contract NAS8-39073 to the Chandra X-ray Center. 
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APPENDIX 



EFFECTS OF VACUUM POLARIZATION 



In § 3.1 and § 3.2 the NS surface emissivity has been computed under the assumption that radiation propagates in 
vacuo outside the star and neglecting the magnetized vacuum birefringence and polarization properties. At the field 
strengths we consider (B < Bqed = m e c 3 /he ~ 4.4 x 10 13 G), this has little affect on the vacuum refractive index 
for which deviations from unity are very small (see below). Yet, radiation propagating in the magnetized vacuum has 
two well-defined polarization states, corresponding to the ordinary (O) and extraordinary (X) modes, even if the two 
modes propagate at very nearly the same speed (no — nx — 1). Therefore, in principle the entire formalism should be 
generalized to account for non-scalar absorption and emission coefficients. In order to quantify this effect, we proceed as 
follows. We maintain the same cartesian frame introduced in § 3.1 and, at fixed magnetic co-latitude 9, we first consider 
incident radiation with polarization mode s (here and in the following s stands for either O or A). The vacuum dielectric 
tensor is expressed as (e.g. Meszaros 1992; Heyl & Hernquist 1997) 



a + q sin a q sin a cos a 

a 

q sin a cos a a + q cos 2 a 

Suitable expressions for a, q have been given by Ho & Lai (2003). In the weak field limit B < Bqed it is 

awl- 25 v , q w 7S V , S v = ^—b 2 v 

457T 

where aj? = 1/137 and by = B/Bqed, while for B > Bqed it is 

OLp 



1 + 



457T 



Olf 

^-45^ 

The two unit propagation eigenmodes are 



1.195- § In by - -!- (0.8553 + In b v ) - 

3 by 2b v 

-\b v + 1.272 - -!- (0.3070 + In fey) - 0.7003-i- 



(Al) 

(A2) 
(A3) 



bv 



ex 
eo 



1 



k x b 



sin 6 sin 5 



(cos a sin i sin (3 , — cos a sin i cos (3 + sin a cos i , — sin a sin i cos (3) 



b — n Q cos Sk 



y/l + n' (n - 2) cos 2 6 ^/l + n (n A - 2) cos 2 6 
—n cos (5 sin i sin f3 , cos a — n Q cos 8 cos i) 



(sin a — n cos 5 sin i cos (3 , 



(A4) 



with cos 5 = k • b = cos i cos a + sin i sin a cos (3. Series expansions for the refractive indices of the two modes have been 
given by Heyl & Hernquist (1997). In the weak-field limit it is 



no 



n x 



1 ap ■ 2 x 
1 — sm o 

Air 

OLE . 2 x 

1 + -— sin o 

Air 



!6 , 2 64 , 



0(b e v ) 



O 



\4irJ 



b v - 0.53 (6B e - 5B 4 ) b v + 0(b v ) 



while the corresponding expressions in the strong-field limit are 



+ o 



/ttfV 
V4tt/ 



(A5) 



n 



n x 



1 , ap ■ 2 
1 + — sm 

47T 

a F . 2 

1 + — sm 2 

47T 



--(ln6 y + l 



hiTr) — + 
by 



3 
+ 



1 2 
8 In A - - - - lE 



ln7r 



+ o 

7T 2 

18 " 



Utt/ 
2 -In b 



/a F y 



(A6) 



1 



here B„ are the Bernoulli numbers, ln^4 ~ 0.2488 and je is the Euler constant. 

At the surface an incident electromagnetic wave, described by its electric field E s , wave vector k s and refraction index 
n s , is partly reflected and partly refracted. Because of the birefringence of both media (the vacuum and the solid star 
crust) this gives rise to two refracted and two reflected waves. The latter are again either X- or O-polarized with electric 
field E" r where, once more, r = O, X 7 . In the following we ignore the small (0\no — nx\) difference in the angle of 
reflection, and assume that the two reflected waves propagate in the same direction fixed by the angles i, [3 + it. 

7 Strictly speaking, the reflected wave is not necessary linearly polarized; however it will separate in the two allowed polarization states after 
propagating a distance l v w 2nc/(uj\no —nx\) ~ 10~ 4 cm for X-ray energies and B « Bqed (see e.g. Chanan, Novick & Silver 1979; Meszaros 
1992). The total intensity of the reflected radiation will be the same far away from the source, despite the change in the polarization state. 
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Since we are not interested in studying the polarization state of the emergent radiation, but only its spectral distribution, 
it is useful to introduce the two quantities 

P "- s = |E.|» (A?) 

which represent the ratios between the reflected and incident intensities when the incident wave is either X- or O- 
polarized (intensities are proportional to squared amplitudes of the electric fields to first approximation 8 ). The absorption 
coefficients corresponding to the two incident modes are a U;S = 1 — and, again, from Kirchhoff's law we get the total 
emissivity j u — a w P w (T), where a w = (a w ,o + &u,x)/2. The monochromatic and total flux are computed again by 
performing the integrals in eqs. (9) and (10). 

Inside the star, each refracted wave splits into an ordinary (E' sl ,\i' sl ) and an extraordinary (E(, 2 ,k' s2 ) mode. In 
order to compute p UtS , we then proceed exactly as discussed in § 3.1 by solving the dispersion relation and computing the 
refractive index n,, i = 1,2, for the two refracted modes. Since n s now explicitly appears in Snell's law (n = n s sini/ sin©), 
we have to solve the dispersion relation twice 

n 4 (P + v sin 2 a) + n 2 (gv - 2PS + u sin 2 a) + PRE + gu = 

n s sinisin(2a) cos/3(n 2 - n 2 sin 2 i) 1 / 2 ^ + n 2 v) , (A8) 

where g = n 2 sin 2 i[l — sin 2 a(l + cos 2 (3)\ and all other quantities are the same as in eq. (15). 

Once the refractive indices are known, we solve the wave equations \s y ij{nm)E' sm: j = 0, where E' sm j are the cartesian 
components of E' s m , obtaining the two ratios E' s ^ mx / E' s ^ mz and E' s ^ my l E' s ^ mz . The resulting expressions are 



E' 
~E' 



— <^m — 



— n m n sin i sin (3 cos [3 — iDn 2 sin i cos a + iDn s cos (3 sin a sin i \ n 2 



n 2 sin 2 i 



— n s sin 

,2 -2 



E' 

- o m - 



in/3 (P — S) sin a cos a sin i \J — n 2 sin 2 i 
+ n 2 sin 2 i sin (3 cos (3 (P cos 2 a + S sin 2 a) + iD cos aP] 

x —n 2 m n s sin i \J — n 2 sin 2 i sin (3 + iDsman 2 n — iDn 2 sin a sin 2 i cos 2 (3 

— iDn s cos a cos (3 sin i \J n 2 ^ — n 2 sin 2 i + n s sin i \J n 2 ^ — sin 2 i [sin 

+ sin (3 sin 2 a (P — S)] — (P — S)n 2 sin a cos a sin 2 i sin /? cos (3 — iD sin aP] 

a m ("s sin 2 * s i n P cos (3 — iD cos a) + n s sin /3 sin i \J n 2 ^ — n 2 sin 2 i + iD sin a 



x (n 2 sin 2 /3 sin 2 i 



ST 1 



(A9) 



(A10) 



The generalization of the Fresnel equations (e.g. Jackson 1975 and B80) to the present case gives 



S E S ± + "r-E^rl 



B m E' 



n s E s j_ - n r E. 



HI E' 



/ j ll r-^ s ,r± — / j w mU m i-J smz 



^H-E<r||=E^< 



B -ll+E<r||=E£^ 



(All) 



In the previous expressions the components of E s , E" r are parallel and orthogonal to the plane of incidence, w n 
\jn 2 m - n 2 sin 2 ij cos i, A m = b m sin j3 - a m cos j3, B m = b m cos (3 + a m sin j3, C m = a m e' 31 + b m e' 32 + e' 33 and e - fe e^ c 



The components of the electric field of the incident and reflected waves can be expressed in terms of the amplitudes E s 
and E" as 



E s ± — f s E s , E s » = g s E s 



P, 



r± 



KK 



p. 



r," J?" 

9rE« r 



(A12) 



where 



fx = cos 7, g x = sin 7, cos 7 



- cos a sin i + cos f3 sin a cos i 
sin S 



8 We neglect deviations between the directions of the two reflected waves and that of the corresponding time-averaged energy fluxes, the latter 
defined by the Poynting vectors. 
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fo = cos£, go = sin£, cos£ 



sin (3 sin a 



^/l + n (n - 2) cos 2 S 



(A13) 



The expressions for f x , /q, g x , q'q are obtained from the previous ones by replacing (3 with j3 + it. 
Inserting (A12) into equations (All) finally gives the amplitude ratios 
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(A15) 



Again, when only one refracted wave (labeled "1" for convenience) survives, the previous calculation yields (f s ,g s ^ 0) 



f s w\ — \ g s Ci cosi — Ai sini 



Ko 



Kx 



n x fx w i + 1 9x cos i + Ai sin i 
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9o 
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9x 
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5o 
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(A16) 



from which p WjS follows. Similar expressions can be derived in the case in which either f s or g s vanishes. 

The absorption coefficients a w s = 1 — p w . s have been computed numerically in the relevant angular ranges following the 
procedure outlined above and the results have been used to evaluate j u and f^ (see eq. [9]). The quantity f^/B^iT) has 
then been compared with its value computed by neglecting vacuum polarization. We repeated the comparison for different 
values of magnetic field; some examples are shown in Figure 10. Fractional corrections turn out to be always negligible, 
being at most ~ 1CP 2 below 2 kcV. Vacuum corrections enter the dielectric tensor via the two quantities a — 1 and q, and, 
even for B p ~ 5x 10 13 G, the largest value we consider in our model, it is \a— 1| ~ q « 10~ 4 . The corresponding deviation 
of the refractive index from unity in the vacuum outside the star is indeed negligible. Therefore, when we compute the 
total reflectivity we are superimposing two incident modes perpendicular to each other. By choosing them in the plane 
perpendicular and parallel to the plane of incidence or, as in this case, with ordinary and extraordinary polarization, is 
not important as long as they travel at nearly the same speed. 



Table 1 

Isolated Neutron Stars Parameters 



Source T bb (eV) a B (10 12 G) b Refs. c 



RX J1856. 5-3754 


61.1 ±0.3 




1, 


2 


RX J0720.4-3125 


86.0 ±0.6 


21.3 ±0.1 


3, 


4 


Vela 


128.4 ±7 


3.3 


5, 


6 


Geminga 


48.3t 6 9 j 


1.5 


7, 


8 


PSR 0656+14 


69.0 ±2.5 


4.7 


9, 


6 


PSR 1055-52 


68.Hi° 7 ; 2 2 


1.1 


10 


,6 



a Errors refer to 2a confidence level 

b As computed from the spin-down formula; period 
derivative is very accurate for the three radiopulsars 
and for Geminga. Errors refer to 90% confidence level 
for RX J0720.4-3125. 

c [1] Burwitz at al. (2001); [2] Drake et al. (2002); 
[3] Pacrcls ct al. (2001); [4] Zane et al. (2002); [5] 
Pavlov et al. (2001); [6] Taylor, Manchester & Lyne 
(1993); [7] Halpern & Wang (1997); [8] Bignami & 
Caraveo (1996); [9] Marshall & Schulz (2002); [10] 
Greiveldinger et al. (1996) 
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12 3 

log B 12 



Fig. 1. — The critical temperature for H and Fe as a function of the magnetic field. Condensation is possible in the shaded region for Fe 
and in the cross-hatched region for H. The full circles with error bars mark the position of five cool, isolated neutron stars (see table 1). The 
horizontal line is drawn in correspondence to the color temperature of RX J1856. 5-3754. 




-1.5 -1.0 -0.5 0.0 

log Energy (keV) 



Fig. 2. — The monochromatic absorption coefficient as a function of energy for B = 10 12 G and different values of the magnetic field angle: 
curves are from top to bottom for 2a/n = 0.05 , 0.2 , 0.4 , 0.6 , 0.8 , 0.95. The plasma frequency given by u> p fi has been used here. 



17 




Log Energy (keV) 

Fig. 3. — The total absorption coefficient averaged over the star surface, F U /B U (T), as a function of energy and for different values of the 
magnetic field: B p = 5 X 10 12 G (full line), 10 13 G (dotted line), 2 X 10 13 G (dashed line), 5 X 10 13 G (dash-dotted line). The two sets of 
curves correspond to models without and with electron-phonon damping accounted for (the latter are evaluated assuming T = 10 6 K). The 
plasma frequency LO Pt o has been used here. 
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Fig. 4. — Ratio of the emitted to the blackbody power in the 0.1—2 keV band for different values of the plasma frequency. Circles refer 
to Bp = 3 X 10 13 G and diamonds to B p = 5 X 10 13 G. Filled and open symbols are for the uniform/meridional variation temperature 
distributions respectively. In the latter case we assumed a profile T 4 (0) = Tf urf [K + (4 - K) cos 2 0]/{[l - 0.47(1 - K)](l + 3 cos 2 9)} with 
K = 10~ 4 (e.g. Greenstein & Hartke 1983; Possenti, Mereghetti & Colpi 1996). 
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Log Energy (keV) 



Fig. 5. — Same as in Figure 2, but with electron-phonon damping accounted for and B = 5x 10 13 G, T = 10 6 K; curves are, from top to 
bottom for 2a/ir = 0.05 , 0.2 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.95. 
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log Energy (keV) log Energy (keV) 



Fig. 6. — Upper panels: Re(n„) with (right) and without (left) electron-phonon damping accounted for. The refractive indices are the 
solutions of the dispersion relation for Kb 43°, a = 0, and B = 5 X 10 12 G (solid line), B = 5 X 10 13 G (dash-dotted line). Lower panels: 
same for the imaginary parts of n m . The dashed line represents our rejection criterion, Im(n) > 0.01 (see text). 
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Log Energy (keV) 



Fig. 7. — Same as in Figure 3 for B p = 5 X 10 G and different values of the rejection limit: Im(n) = 0.01 (solid line), Im(n) = 0.05 
(dashed line) and Im(n) = 0.005 (dotted line). The dash-dotted line shows the surface emissivity for the meridional temperature distribution 
(Im(n) = 0.01; see the caption of Figure 4) and is to be compared with the solid curve. 




log Energy (keV) log Energy (keV) 

Fig. 8. — The emitted spectrum in the cold plasma limit for B v = 2 X 10 13 G and T aur f = 10 6 K. Left panel: uniform surface temperature; 
right panel: meridional temperature variation as defined in the caption of Figure 4. The dashed line is the blackbody at T sur f and the 
dash-dotted line the blackbody which best-fits the calculated spectrum in the 0.1-2 keV range. The two models shown in each panel are 
computed for u> p = ui p fi (upper solid curve) and 2.45oj Pj o (lower solid curve). Spectra are at the star surface and no red-shift correction has 
been applied. 
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Fig. 10. — Fractional difference in the computed values of / B^(T) with and without vacuum polarization taken into account (see 
Appendix A) at different energies. Solid and dotted lines arc for B p = 3 and 6 X 10 13 G respectively. 



